Cluster Dynamics for Randomly Frustrated Systems with Finite Connectivity 
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In simulations of some infinite range spin glass systems with finite connectivity, it is found that 
for any resonable computational time, the saturated energy per spin that is achieved by a cluster 
algorithm is lowered in comparison to that achieved by Metropolis dynamics. The gap between 
the average energies obtained from these two dynamics is robust with respect to variations of the 
annealing schedule. For some probability distribution of the interactions the ground state energy is 
calculated analytically within the replica symmetry assumption and is found to be saturated by a 
cluster algorithm. 



Many systems composed of a macroscopic number of 
interacting elements share the property of being compu- 
tationally difficult. By this we mean that their relaxation 
time, as well as the time scales related to the system in- 
vestigation, grow very fast with the size of the system [Q . 
This concerns both the actual dynamics of the physical 
systems as well as pseudo dynamics used in computer 
simulation Q|. Such systems appear in many physical 
fields, from statistical mechanics, to the study of quan- 
tum field theories (for review sec 

A typical example of such difficulties is the critical 
slowing down at second order phase transitions. This 
phenomenon is simply the divergence of the relaxation 
time as the critical point is approached. Consequently, 
the typical time needed to produce a large Boltzmann set 
of decorrelated configurations diverges and the standard 
local Monte Carlo simulation methods become inefficient. 

To solve such problems, multiscale-cluster-type algo- 
rithms were devised, and the entire subject of global 
collective dynamics attracted considerable attention. It 
is generally believed j|] that for several classes of sys- 
tems, multiscale methods may overcome the slowing- 
down problems. Moreover, the multiscale algorithms are 
conceptually important insofar as they encode the un- 
derstanding of the relevant large-scale physics. In par- 
ticular, these procedures isolate the relevant degrees of 
freedom and act directly on them, in a manner con- 
sistent both with their effective macroscopic dynamics, 
and with the basic interactions that define the system 
at the microscopic scale (for a review see Ref. [^). The 
cluster-multiscale methods have been applied successfully 
to many fields in physics (second-order phase transitions, 
disordered systems, quantum field theories, fermions in a 
gauge background, quantum gravity, and more and 
in many cases a dramatic acceleration of the numerical 
simulation was achieved. 

However, the general applicability of the multi scale 
methods is in question. The situation is particularly un- 
clear for models with complex energy landscape. Their 
physical properties are notoriously hard to investigate, 
especially in the low-temperature phase. Some of the 
most important families of systems presenting such dif- 
ficulties are the randomly frustrated systems (RFS) and 
in particular spin-glass (SG) systems. 

The study of SG has attracted wide research activ- 



ity over the last two decades (for review see yj5|) and 
their theoretical understanding goes beyond its original 
scope of understanding experimental results of real phys- 
ical systems . The progress of the statistical mechanics 
methods related to SG contributed to the understanding 
of a wide variety of other disordered systems. One of 
the most promising directions is to apply the SG knowl- 
edge to the study of hard optimization problems belong- 
ing to the NP class This relationship goes beyond 
a mere analogy, and the task of determining the optima 
of a problem can be rigorously mapped in to finding the 
ground state (GS) of the analogous SG system. A typical 
SG system presents a complex energy landscape consist- 
ing of many local minima, separated by huge barriers 
that scale with the size of the system, and lead to an 
infinite hierarchy of exponentially divergent relaxation 
time scales Therefore, besides the problem of proper 
sampling in simulations at low temperatures, the system 
has a tendency to get stuck in local minima that pre- 
vent the measurement of equilibrium properties within a 
reasonable time. Furthermore, the simulated annealing 
technique that prevents some systems from getting stuck 
in local minima, failed to provide a complete solution in 
SG systems. 

On the other hand, the general question of the exis- 
tence of efficient cluster algorithms (CA's) for RFS is still 
an open problem even though a number of successes have 
been achieved for some particular systems. For instance, 
Kandel et al. found an efficient CA for a special case of 
a fully frustrated system on a square lattice, but with the 
lack of randomness. In the special case of two-dimensions 
Ising SG, in which its GS can be found in a polynomial 
time in contrast to the NP feature of the general SG, 
Swendsen and Wang Q developed their "replica" algo- 
rithm. Nevertheless, in general, the cluster algorithms 
are unable to identify the important large-scale degrees 
of freedom in RFS, and therefore they show no improve- 
ment on simple local algorithms. 

Hence, an efficient CA for RFS, if any, will enhance 
our understanding of the low-temperature physics of SG 
systems, and its consequences on related problems, such 
as efficient heuristics algorithms for solving NP problems. 
In particular, finite connectivity models at low tempera- 
tures are directly connected to graph partitioning . 
This is the problem of dividing a given graph into sub- 
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graphs, with minimum connections between them. Be- 
side these apphcations, an efficient CA for SG may serve 
as a tool for the understanding of rephca symmetry brak- 
ing (RSB) and its properties. Note that a quantitative 
fitting to Parisi's picture, and the existence of RSB in fi- 
nite connectivity models, are still open and controversial 
questions. 

In this letter we present a step towards understanding 
the applicability and the limitation of CA for RFS. We 
consider an Ising system described by the Hamiltonian 



(1) 



where Si — H (i — 1, 
bution of the links is 



■ ,N) and the probability distri- 



F(J,,) = (1 - c/N)6{J,,) + {c/N)f{J,,). 



(2) 



This model is known as a highly diluted system with 
finite connectivity, since the probability for a spin 
with connectivity k follows the Poisson distribution: 

exp(— c)/fc! with average connectivity c, which is taken 
to be 0(1), and f{Jij) is the distribution of the sur- 
viving links after the dilution. In the present work we 
would consider the following types of unbiased [/( ) = 
/(— Jy)] distribution: (a) Gaussian distribution, (b) 
Jij = ±1, and (c) a special case in which the links get 4 
values: Jij = ±1, Jij = ±e. 

This model with Jij = ±1 was systematically stud- 
ied near the glass transition temperature by Viana and 
Bary and at low temperatures by Kanter and Som- 
polinsky ||l^ and by Mezard and Parisi [0. The self- 
consistent description of the low temperatures is based 
on the probability distribution of the local field defined 
by hi = Ttanh"^ < 5'^ >t- Physically, this field is the 
first excitation, namely, in the limit T — > 0, \hi\ is the 
minimum energy cost for flipping the ith spin from its GS 
by the "best" reorganization of the system. This local 
field is in truth an oxymoron, since it depends on global 
properties of the cluster, the exchange field '^j Jijfnj, 
on the other hand, is truly a local property depending 
on the local connectivity (note that \hi\ < | Jijirijl). 
In a simple ferromagnetic case (that is J=l) as T — > 
the distribution of the local fields reduces to a discrete 
spectrum ||lo|,[ll| 



Pi = (cP)'exp(-cP)//! / = 0, ...,oo. 
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where P is the fraction of the spins belonging to the 
macroscopic cluster. The quantities Pi, characterize 
global properties, for instance, in the ferromagnetic case 
Pq is the fraction of spins belonging to the finite clusters 
and Pi is the fraction of spins that can be disconnected 
from the infinite cluster by cutting only one link. 

Many geometric properties of this system are well un- 
derstood |lO|JTl| ]. In particular, the system undergoes 
a percolation transition at c = 1. The maximal clus- 
ter is of 0{logN) for c < 1, 0(A^2/3) at c = 1, and 



of 0{N) for c > 1 where its size is explicitly given by 
P = 1 - Po = 1 - exp(-cP). 

The topological structure of diluted models makes 
them good candidates for employing CA's. CA's usually 
consist of two main steps. First, blocks are constructed 
stochastically from many single elements where the links 
between them have been " frozen" according to their ten- 
dency to act coherently, and the other links deleted. Sec- 
ond, updates are performed in which entire blocks flip 
rigidly, in such a way that the Gibbs distribution is still 
fulfilled. Accordingly, the lattice splits into a set of clus- 
ters each formed by sites that can be linked by a chain of 
frozen links. For general SG systems this procedure fails 
because it freezes the entire lattice into a single block. 
However, for a highly dilute lattice the additional dele- 
tions may be just enough to actually split the lattice into 
disjointed blocks. In this way, the CA are effective in lo- 
cating and acting upon large regions that interact weakly 
with the rest of the configuration. In addition, the CA 
may have an efficient way of pinning the frustration to 
weak links. The global decisions made by the CA on 
large-scale degrees of freedom are directly connected to 
the structure of the field that presents global features. 

Simulations on the model, Eqs. (1) and (2), and with 
a Gaussian distribution for the links, /(J) oc exp(— J^), 
were carried out using both local (Metropolis) and clus- 
ter (Wolff ||l^) dynamics. The simulations were carried 
out for various connectivity values, and at temperatures 
below the glassy transition Tc- The size of the system 
was between 1000 and 5000, and the results were aver- 
aged over at least ten different samples. A typical result 
is presented in Fig. 1, and for all times monitored (up to 
100 000 MCS per spin), one can clearly see a gap between 
the energies of the two dynamics. Clearly, in the limit 
of an extremely long time, this gap should vanish, but 
in practice it can be considered as two different energy 
levels. 
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FIG. 1. E{t) per spin for the Gaussian distribution with 
average | J| = 1, c = 2, = 5000, and T ~ O.ST^. The solid 
and the dotted hues indicate Wolff and Metropolis dynamics, 
respectively. The time scale is in MC steps per spin. 
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An interesting question is whether this gap is robust 
even under an annealed schedule in which the temper- 
ature is gradually lowered. This is done to avoid be- 
ing caught in " false minima" , which is typical for low- 
temperature simulations, and also enable the CA to act 
on clusters at various scales. Indeed, the results were im- 
proved for both dynamics, but the gap still exists [ p^ . 
On one hand it is clear that in practice at low temper- 
ature, the energy Ec reached by the CA is lower than 
that of the local algorithm Ei. However, on the other 
hand, this improvement should be evaluated with re- 
spect to the true equilibrium energy, E. Quantitatively, 
ii \Ec — Ei\ <C \E — Ei\, then even after the improvement 
by the CA, the system is far from equilibrium. Hence, the 
calculation of the true GS is important in order to evalu- 
ate the improvement achieved by the CA. Unfortunately, 
the analytical calculation of the GS energy for the Gaus- 
sian distribution appears to be a very difficult task since 
the local field is a continuous variable . We therefore 
need a model whose GS energy can be calculated analyt- 
ically and where the difference between Metropolis and 
CA is enhanced. We find these two features in what we 
shall call the Weak-Strong (WS) model. In this model 
some of the links are much stronger (in their absolute 
value) than others, and explicitly the link distribution is 
given by 

fiJ) = l[6iJ-e) + S{J + e)] 



(1-a) 



[6{J ^ 1) + 6iJ + 1)]. 



In the first set of runs we chose the connectivity c and 
the fraction of the strong links (1 — a), such that the 
density of the strong links by themselves is below the 
percolation threshold, eg = c(l — a) < 1. It is clear that 
in this situation all the strong links are unfrustrated, and 
the frustration is located only on the weak links. This 
framework simplifies the analytical treatment and only 
the energy of the weak hnks, Ew, is considered. In Fig. 
2 one can see (up to our running time) a large steady 
difference (35%) in the energy between the local dynam- 
ics and that of the cluster dynamics. This difference is 
due to the fact that the local dynamics is totally stuck, 
since the probability of flipping a cluster consisting of 
strong links is practically zero for local algorithms. On 
the other hand, the CA "knows" how to deal with the 
strong link structures by considering them as only "one 
degree of freedom" for each cluster. In other words, the 
CA is extremely efficient in the following problem: "How 
to arrange the weak links in the environment of strong 
links" . Note that the origin of this effect is the same as 
in the Gaussian case. In the WS case, however, it is am- 
plified, due to the special discrete link distribution. In 
order to emphasize the CA features, we performed mea- 
surements using simulated annealing for both dynamics. 
This was performed over a wide range of temperatures, 
enabling us to first arrange the strong links, and then 
lowered the temperature to the scale of the weak links. 
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FIG. 2. Ew{t) per spin for the WS case, c = 2, a = 0.7, 
— 5000, T = 0.5 Jw, the strong and the weak links was 
scaled to 100 and 1, respectively. The solid and the dotted 
lines indicate Wolff and Metropolis dynamics, respectively. 
Inset; The first 50 steps. 

Figure 3 demonstrates one of the examined schedules 
from which one can conclude that the "gap" (40%) is 
robust (or even increases) with respect to the annealing 
schedule. This provides strong evidence that the CA in- 
deed do overcome the local difficulties. In order to put 
the CA improvement in the overall picture of the true 
equilibrium, one must calculate the GS energy. After 
some calculation, one can show that the GS energy of 
the weak links within the replica symmetry assumption 
is given by 
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Ew = -^caP^e + -c(l - a) 5](1 - 4a^)e 



k=Q 



(4) 



where = ^ + ZLi ^ = Y.h=-oo \h\P{h). Note 
that the energy of the strong links, Es = — |c(l — a), is 
eliminated from Eq. (^). 

The explicit value of Ew depends on the local field, 
P{h), which in general is difficult to calculate. Never- 
theless, in the limit eg < 1 and e <C -^j which implies 
h ^ 1 (there are no "order 1" local fields), the following 
distribution of local fields can be assumed 



Pih) = {1-Q)dih) 



PiS{h-le), 



(5) 



where P/ ~ P^i , and Q = 1 — Pq the fraction of the frozen 
(zero entropy at T = 0) spins. The quantitative form of 
P{h) can be determined from the following self-consistent 
equation: 



P{h) = e 



dy 



1=1 



iyle 



-iyle 



(6) 



In principle, one must solve infinitely many coupled non- 
linear equations, obtained from the expansion of Eq. m). 
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However, assuming that Pi decays with I, we assume 
Pi = ioT I > I', and then from the comparison of the 
expansion of Eq. (||) and Eq. ^ one can derive I' nonhn- 
ear equations. For the range of c and a in our simulations 
we take /' = 8, which gives a negligible error. For com- 
parison we calculated Pi in the J — ±1 case in which we 
obtain 



Pi = exp(-c(3)/|/|(cQ), 



(7) 



where Ii{x) is the modified Bessel function. A graph of 
Pi for the two cases is presented in the inset of Fig. 3. 

One should note that after scaling e to 1, the Pi for the 
two cases is very close. However, the exchange field is a 
different story. In the J = ±1 case the exchange field 
of a spin is usually not far from its local field (around 
the number of its neighbors), but for the WS case the 
local field is 0(e) while the exchange field may be 0{1). 
The fact that Pi is much smaller leads us to hope that 
the dynamics that prefer a global "decision" are indeed 
superior to the local one. 

The type of the mean-field solution we derived is known 
to be unstable |Q. However, in Fig. 3 one can see that 
the analytical GS obtained from Eqs. ^ and (H), is in 
very good agreement with the averaged GS energy ob- 
tained by the CA. 
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FIG. 3. Ewit) per spin for the WS case, c = 2 and a = 0.7. 
The annealing schedule range is T G 123 — 0, with AT = 10 
for T > 3 and AT = 0.1 for T < 3. The jump in Ew at 
T = 3 is observed since we are getting close to the scale of 
the weak link and NEw become extensive. The solid and 
the dotted lines indicate Wolff and Metropolis dynamics, re- 
spectively. The horizontal line denotes the analytical GS in 
Eq. (4). Inset: Pi, for the WS case (solid), and J = ±1 case 
(dotted). 

The scope of our analytical calculation is for the 
cs < 1 case (below the percolation threshold for the 
strong links). However, we also performed simulations 
for Cs > 1, and the results for c G (2,5) show a simi- 



lar picture. Thus, the superiority of CA over the local 
dynamics is shown clearly. However, apriori, one could 
gain the impression that this was achieved only by some 
simple geometrical reduction of the system. In order to 
check this point we performed another type of simula- 
tions (which will be reported in detail elsewhere [|l4|) in 
which the algorithm explicitly reduced all the trees, the 
linear chains, and the self-loops. This was carried out 
in a manner that is proved to be energetically optimal. 
Nevertheless, we show that the CA advantage goes be- 
yond these simple reductions. This result completes the 
picture of CA superiority, within the scope of the present 
letter. There are many questions that still remain open. 
First, what happens at high connectivity? It is clear that 
our CA superiority decreases as c grows. One can investi- 
gate this point further, and classify the relevant windows 
of parameters. Moreover, one can ask if it is possible to 
overcome this limitation by a new type of CA, namely, a 
CA that goes beyond the ability to act on blocks having a 
considerably weak interaction with the rest of the system. 
Second, an interesting question is whether CA can enter 
the RSB region. Our results do not clarify this point, and 
indeed our results are consistent with the RS GS energy. 
However, the existence of RSB in the examined systems 
is still in question. 
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